Lab 9: Spaceborne Lidar
Anthony Stewart
2023-05-22
Objectives:
Introduction to full waveform lidar
Introduction to Global Ecosystem Dynamics Investigation (GEDI) Lidar
Introduction to The Ice, Cloud and land Elevation Satellite-2, or ICESat-2
Downloading and using GEDI data
Characterizing forested landscapes with full waveform lidar
Data and Software:
LAB9_GEDI data
rGEDI: An R Package for NASA’s Global Ecosystem Dynamics Investigation (GEDI) Data Visualization and Processing https://cran.r-project.org/web/packages/rGEDI/vignettes/tutorial.htmlArcGIS Pro
What you will turn in:
Welcome to Lab 9 for ESRM433/SEFS533 This lab is all about Spaceborn lidar. We have covered aerial lidar, and terrestrial lidar, now we will talk about spaceborn lidar scanning (SLS). Two prominent spaceborn lidar sensors are ICESat and GEDI.
PART 1a GEDI background
For this lab, we are going to focus on GEDI data. The GEDI sensor was launched to the International Space Station (ISS) in December of 2018 and is slated for a two year mission to scan terrestrial ecosystems. The data collected is limited to the orbit of the ISS so data isn’t available for latitudes greater than 51.6o.
For the following questions refer to: - GEDI’s website https://gedi.umd.edu/
GEDI videos - https://www.youtube.com/watch?v=SSdDPFfUVIo - https://www.youtube.com/watch?v=wxgrxvAKpTo - https://www.youtube.com/watch?v=XFlm-TmhvjM
Any other online resources you can find for GEDI and full waveform lidar
The lidar data from GEDI retains it’s full waveform. ALS data is also based on the amount of energy returned from a laser pulse, but the full waveform is discarded and only the peaks in the waveform are preserved. These peaks are what we call returns when using terms like “discreate return lidar” vs “full wavefrom lidar”
The beam divergence of lasers is extremely important to keep in mind when discussion SLS. Beam divergence is in direct relationship to the distance the laser pulse travels. For TLS, objects are much closer to the scanner so beam divergence is likely to be a few mm to a few cm depending on the distance. For ALS, beam divergence can be ~10cm to ~30cm. This divergence is what allows for multiple returns from a single pulse.
QUESTION 1: Is GEDI data free to the public? Where can you get GEDI data?
QUESTION 2: What is the footprint (resolution) of GEDI laser pulses? This is the same as beam divergence.
QUESITON 3: How many beams of data is GEDI collecting? What is their spacing?
QUESTION 4: Is every square meter of the earths surface getting sampled by GEDI?
QUESTION 5: What are the three science questions for GEDI?
QUESTION 6: For Data Products, there are multiple Algorithm Theoretical Basis Documents (ATBDs) for GEDI data. We are most interested in L1B, L2A, & L2B. What is the name of the data product for those 3 ATBDs?
PART 1b ICESat background
ICESat is like GEDI in many respects. An orbiting lidar sensor in space. A major difference is that GEDI is on the international space station while ICESat has its own dedicated satellite in a polar orbit. Take a moment to review the information here about ICESat: https://www.nasa.gov/content/goddard/about-icesat-2
QUESTION 7: What are some similarities and some differences between GEDI and ICESat2?
QUESITON 8: What wavelength and color is the ICESat2 laser? What wavelength and color is the GEDI laser?
QUESTION 9: What is the spacing of the ICESat2 laser pulses?
PART 2: Getting and processing GEDI data
Most of the following tutorial was taken from: https://cran.r-project.org/web/packages/rGEDI/vignettes/tutorial.html
An important part of this class is to demonstrate how accessible lidar data is for researchers. This is certainly true about GEDI data, however, GEDI data downloads are huge. The smallest bundles are about 6GB. This is in the highly compressed .h5 format. H5 is a Hierarchical Data format and if you want more info check out: https://www.neonscience.org/about-hdf5
Because of the size constraints, you will be provided with a small snippet of GEDI data to work with. A package for R has been developed to facilitate working with GEDI data, rGEDI.
rGEDI WILL NOT WORK ON MacOS! The package hasn’t seen many updates and will only work on Windows. Therefore you may need to work through Madrona/SEFS Spatial
library(devtools)
devtools::install_git("https://github.com/carlos-alberto-silva/rGEDI", dependencies = T)
This may take a while to run, especially if you say yes to updating
all the packages. I would suggest saying no or
3 for now unless you run into a problem related to package
versions.
library(rGEDI)
library(sf)
library(terra)study_area <- st_read("Lab9/Lab9Data/study_area.shp")
study_extent <- st_bbox(study_area)ul_lat<- study_extent["ymax"]
lr_lat<- study_extent["ymin"]
ul_lon<- study_extent["xmin"]
lr_lon<- study_extent["xmax"]daterange=c("2021-06-01","2021-08-01")?gedifinder
gLevel1B<-gedifinder(product="GEDI01_B",ul_lat, ul_lon, lr_lat, lr_lon,version="002",daterange=daterange)
gLevel2A<-gedifinder(product="GEDI02_A",ul_lat, ul_lon, lr_lat, lr_lon,version="002",daterange=daterange)
gLevel2B<-gedifinder(product="GEDI02_B",ul_lat, ul_lon, lr_lat, lr_lon,version="002",daterange=daterange)gLevel1B
gLevel2A
gLevel2Bset an output directory for gedi
outdir <- ("Lab9/Lab9Data/GEDI_DL/")Need a NASA EarthData Explorer account
gediDownload(filepath=gLevel1B,outdir=outdir)
gediDownload(filepath=gLevel2A,outdir=outdir)
gediDownload(filepath=gLevel2B,outdir=outdir)Just doing the file reading
You now have GEDI data loaded into R! There is a lot of data packed into those h5 files. The first thing that we should do is to get the location information for the pulses.
gedilevel1b <- readLevel1B("Lab9/Lab9Data/LAB9_GEDI/PF_level1b_clip.h5")
gedilevel2a <- readLevel2A("Lab9/Lab9Data/LAB9_GEDI/PF_level2a_clip.h5")
gedilevel2b <- readLevel2B("Lab9/Lab9Data/LAB9_GEDI/PF_level2b_clip.h5")That pulls the lat and long for each shot (pulse) from the h5 file. The header of the new file shows you a snippet of the data. You have shot_number (every pulse has it’s own ID number) and the lat, long, and elevation.
?getLevel1BGeo## starting httpd help server ... donelevel1BGeo <- getLevel1BGeo(level1b = gedilevel1b, select = c("elevation_bin0"))##
|
| | 0%
|
|= | 2%
|
|=== | 4%
|
|==== | 6%
|
|====== | 8%
|
|======= | 10%
|
|========= | 12%
|
|========== | 15%
|
|============ | 17%
|
|============= | 19%
|
|=============== | 21%
|
|================ | 23%
|
|================== | 25%
|
|=================== | 27%
|
|==================== | 29%
|
|====================== | 31%
|
|======================= | 33%
|
|========================= | 35%
|
|========================== | 38%
|
|============================ | 40%
|
|============================= | 42%
|
|=============================== | 44%
|
|================================ | 46%
|
|================================== | 48%
|
|=================================== | 50%
|
|==================================== | 52%
|
|====================================== | 54%
|
|======================================= | 56%
|
|========================================= | 58%
|
|========================================== | 60%
|
|============================================ | 62%
|
|============================================= | 65%
|
|=============================================== | 67%
|
|================================================ | 69%
|
|================================================== | 71%
|
|=================================================== | 73%
|
|==================================================== | 75%
|
|====================================================== | 77%
|
|======================================================= | 79%
|
|========================================================= | 81%
|
|========================================================== | 83%
|
|============================================================ | 85%
|
|============================================================= | 88%
|
|=============================================================== | 90%
|
|================================================================ | 92%
|
|================================================================== | 94%
|
|=================================================================== | 96%
|
|===================================================================== | 98%
|
|======================================================================| 100%head(level1BGeo)## shot_number latitude_bin0 latitude_lastbin longitude_bin0
## 1: 72230000300104806 46.89253 46.89253 -122.3994
## 2: 72230000300104807 46.89230 46.89230 -122.3987
## 3: 72230000300104808 46.89208 46.89207 -122.3980
## 4: 72230000300104809 46.89185 46.89184 -122.3974
## 5: 72230000300104810 46.89162 46.89162 -122.3967
## 6: 72230000300104811 46.89140 46.89139 -122.3960
## longitude_lastbin elevation_bin0
## 1: -122.3993 202.2359
## 2: -122.3987 190.0121
## 3: -122.3980 204.2407
## 4: -122.3973 220.8789
## 5: -122.3967 230.3352
## 6: -122.3960 212.5780There are a few steps we need to do to convert the data into a format
that we can output as a shapefile: - Convert shot_number
from ‘interger64’ to ‘character’
level1BGeo$shot_number <- as.character(level1BGeo$shot_number)Convert level1BGeo to a spatial object using the
sf package and giving it a coordinate reference system
(crs) with the EPSG code of 4326
level1BGeo_sf <- st_as_sf(level1BGeo, coords = c("longitude_bin0", "latitude_bin0"), crs = 4326)Now we can finally look at the locations of our individual shots
mapview(level1BGeo_sf, map.types = "Esri.WorldImagery")And we can write the locations as a shapefile (.shp) to
our computer and view it in another program like ArcGIS Pro
sf::write_sf(level1BGeo_sf, "Lab9/Lab9Data/GEDI_level1BGeo.shp")ArcGIS
You should now have a shape file called level1b.shp Open ArcGIS and make a new map project called LAB9. Insert a basemap that is world imagery.
Drag and drop your level1b.shp file into ArcGIS. You know that the diameter of each shot is 25 to 30m. Let’s create a buffer around each point so we get a good idea of the area actually sampled by each shot.
Use the toolbox to search for buffer
Remember, your point size will change depending on your level of zoom. You need the buffers drawn to see the footprint of the shots.
QUESTION 12: Identify 3 different shots:
72230000300104862 - barren or cleared shot
• 72230300300100204 - heavily forested shot
• 72230300300100165 - light to moderate vegetated shot
Include a screenshot of each shot (use the buffer, not the point), AND the coordinates in NAD83/Washington State Plane South of the plot center.
rGEDI
We are now going to compare the full waveform lidar of each shot.
?getLevel1BWF
wfB <- getLevel1BWF(gedilevel1b, shot_number = "72230000300104862")
plot(wfB, polygon = TRUE, type = "l", lwd = 2, col = "forestgreen", xlab= "Waveform Amplitude", ylab = "Elevation(m)", main = "Barren")
QUESTION 13: Include screenshots of your three waveforms. Make sure they are labeled “Barren”, “Forested”, and “Light Vegetation”. Discuss the similarities and differences between the waveforms. Remember, this is the full waveform and is not yet cropped with the ground and tree height identified.
QUESTION 14: You have now seen the output from L1B Give a more thorough description of what exactly L1B is
Lets take a look at the data contained in our L2A sub.
?getLevel2AM
level2AM <- getLevel2AM(gedilevel2a)##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%head(level2AM[, c("beam", "shot_number", "elev_highestreturn", "elev_lowestmode", "rh100")])## beam shot_number elev_highestreturn elev_lowestmode rh100
## 1: BEAM0000 72230000300104806 157.7767 121.6699 36.10
## 2: BEAM0000 72230000300104807 145.7397 122.0679 23.67
## 3: BEAM0000 72230000300104808 159.7441 122.7010 37.04
## 4: BEAM0000 72230000300104809 176.2708 131.9998 44.27
## 5: BEAM0000 72230000300104810 186.4405 155.6543 30.78
## 6: BEAM0000 72230000300104811 164.8258 160.5187 4.30Converting shot_number as “integer64” to “character”
level2AM$shot_number<-as.character(level2AM$shot_number)Converting Elevation and Height Metrics as data.table to
sf spatial object
level2AM_sf<-st_as_sf(level2AM, coords = c("lon_lowestmode", "lat_lowestmode"), crs = 4326)You can export all of the point location and information into an ESRI shape file, just like you did with the L1B data. This time when you click on each point, you will have a list of all the elevation and height metrics.
Export Elevation and Height metrics as an ESRI shapefile
sf::write_sf(level2AM_sf, "Lab9/Lab9Data/GEDI_level2AM.shp")Lets plot the waveform L2A metrics:
?plotWFMetrics
plotWFMetrics(gedilevel1b, gedilevel2a, "72230000300104862", rh=c(25, 50, 75, 90), main = "Barren")##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%
plotWFMetrics(gedilevel1b, gedilevel2a, "72230300300100204", rh=c(25, 50, 75, 90), main = "Forest")##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%
plotWFMetrics(gedilevel1b, gedilevel2a, "72230300300100165", rh=c(25, 50, 75, 90), main = "Light Veg")##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%
For information about specific the algorithms for Canopy Cover and Height refer to: https://lpdaac.usgs.gov/documents/588/GEDI_FCCVPM_ATBD_v1.0.pdf
QUESTION 15: What exactly is the RH metric measuring? Check out: https://gedi.umd.edu/mission/technology/
QUESTION 16: Create plots for your Barren, Forested, and Light vegetation shots. Compare and contrast the waveforms between the different land cover. If there isn’t a clear difference between your plots, sample a new location. Include a screenshot of the three plots.
The last clip of data that you have is L2B. Just as with the other data clips we will bring the data into R and then convert it into a format that can be output as a shapefile:
?getLevel2BVPMQUESTION 17: .
level2BVPM<-getLevel2BVPM(gedilevel2b)##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%head(level2BVPM[,c("beam","shot_number","pai","fhd_normal","omega","pgap_theta","cover")])## beam shot_number pai fhd_normal omega pgap_theta cover
## 1: 0 72230000300104806 2.00382900 3.279090 1 0.3670606 0.63274115
## 2: 0 72230000300104807 0.85066009 3.127365 1 0.6534669 0.34642449
## 3: 0 72230000300104808 2.84563279 3.293170 1 0.2409268 0.75883543
## 4: 0 72230000300104809 4.03077650 3.169244 1 0.1331845 0.86654395
## 5: 0 72230000300104810 2.19123483 3.193733 1 0.3342184 0.66557312
## 6: 0 72230000300104811 0.05784217 1.555367 1 0.9714843 0.02850675Converting shot_number as “integer64” to “character”
level2BVPM$shot_number <- as.character(level2BVPM$shot_number)Converting GEDI Vegetation Profile Biophysical Variables from data.table to sf spatial object
level2BVPM_sf <- st_as_sf(level2BVPM, coords = c("longitude_lastbin", "latitude_lastbin"), crs = 4326)Exporting to ESRI shapefile
write_sf(level2BVPM_sf, "Lab9/Lab9Data/GEDI_level2BVPM.shp")PAI
rGEDI uses a metric called Plant Area Index
(PAI). This is similar to leaf area index but it doesn’t
distinguish between live green leaves and woody tree stems. https://en.wikipedia.org/wiki/Leaf_area_index https://lpdaac.usgs.gov/documents/588/GEDI_FCCVPM_ATBD_v1.0.pdf
PAI differs from height metrics as it relies on the amount of energy that filters through the different layers of forest canopy. A forest with shorter trees but the trees have very full crowns and there is a full midstory layer of vegetation could have a much higher PAI than a forest with very tall trees, but a small crown ratio and little to no midstory vegetation.
Tang, Hao, et al. “Retrieval of vertical LAI profiles over tropical rain forests using waveform lidar at La Selva, Costa Rica.” Remote Sensing of Environment 124 (2012): 242-250.
Zhao, Feng, et al. “Measuring effective leaf area index, foliage profile, and stand height in New England forest stands using a full-waveform ground-based lidar.” Remote Sensing of Environment 115.11 (2011): 2954-2964
LAI and PAI can be expressed as ratios with a value of 1 indicating that the area of plant matter surfaces is equal to the spatial planar area of the sample location. A value less than 1 will have less vegetation, while a value greater than 1 will have more vegetation. There is no upper limit to the lai values but a value approaching 10 would indicate a very dense forest. A value of 0 would be bare ground.
You can look up your PAI values in the biophysical data by sorting
through the table. The easiest way to do it in R is to quickly load in
the library dplyr. This isn’t a remote sensing package,
just a package for working with R data
library(dplyr)
Barren <- level2BVPM |> filter(shot_number == "72230000300104862")
Barren$pai
Forest <- level2BVPM |> filter(shot_number == "72230300300100204")
Forest$pai
Lightveg <- level2BVPM |> filter(shot_number == "72230300300100165")
Lightveg$paiQUESTION 17: What is the PAI for your Barren, Forested, and Light Vegetation shots?
We can create profiles of the PAI for the individual GEDI beams. Think of this as a slice that is taken of the vegetation along the beams path. You can get information about the vegetation height as well as the density of the vegetation. This is the core of calculating the biomass of an area.
let’s get the PAI profiles from the GEDI Level2B
?getLevel2BPAIProfile
level2BPAI_Profile <- getLevel2BPAIProfile(gedilevel2b)##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%head(level2BPAI_Profile[,c("beam", "shot_number","pai_z0_5m","pai_z5_10m")])## beam shot_number pai_z0_5m pai_z5_10m
## 1: BEAM0000 72230000300104806 2.00382900 1.6610067
## 2: BEAM0000 72230000300104807 0.85066009 0.7201815
## 3: BEAM0000 72230000300104808 2.84563279 2.6006696
## 4: BEAM0000 72230000300104809 4.03077650 3.0375919
## 5: BEAM0000 72230000300104810 2.19123483 1.9147891
## 6: BEAM0000 72230000300104811 0.05784217 0.0000000You can create a profile of the different beams from GEDI. You can
determine the name of the beam by selecting a point in ArcGIS just like
how you determined the shot number. You can “normalize” the profile by
removing the elevation from the profile (elev=FALSE).
?plotPAIProfile
# The full power beams are 0101, 0110, 1000, 1011
gPAIprofile <- plotPAIProfile(level2BPAI_Profile, beam = "BEAM0101", elev = FALSE)
gPAIprofile <- plotPAIProfile(level2BPAI_Profile, beam = "BEAM0101", elev = TRUE)
#coverage beams are:0000, 0001, 0010, 0011
gPAIprofile <- plotPAIProfile(level2BPAI_Profile, beam = "BEAM0001", elev = FALSE)
gPAIprofile <- plotPAIProfile(level2BPAI_Profile, beam = "BEAM0001", elev = TRUE)
QUESTION 18: Select a different beam. Include a screen shot of the beam with imagery base map and identify the beam in the image. Include screenshots of the PAI profile with and without the elevation. Caption the images describing what PAI is.
Grid Metrics
Lastly we are going to output a raster of the GEDI data. This is very similar to the grid metrics we used with the ALS data. Consider that we only have spot data for a relatively small proportion of our area. We have to define a grid step size that will each grid cell will include at a minimum one of the GEDI shots. Within each cell, the values of each shot within will used to determine the value assigned to the entire cell. Typical metrics are Max, Min, Mean, and standard deviation.
We need to create a function within R for those metrics:
mySetofMetrics <- function(x){
metrics = list(
min = min(x),
max = max(x),
mean = mean(x), #average
sd = sd(x) #standard deviation
)
}You defined your own metrics using lidR and this is no different.
We are going to grid our metrics into a raster. Make sure to check out ?gridStatsLevel2am. The output is a raster brick like we used before. The resolution is much trickier. The resolution output is in decimal degrees. This will output a different resolution raster depending on the latitude of the sample area. You can use a latitude and longitude calculator to figure out the length of a degree of longitude and latitude at different latitudes. http://www.csgnetwork.com/degreelenllavcalc.html
At our location (latitude ~46.8), 1 degree longitude is ~ 76339m and 1 degree latitude is ~ 111167m.
?gridStatsLevel2AM
rh100metrics <- gridStatsLevel2AM(level2AM = level2AM, func = mySetofMetrics(rh100), res = 0.005)** QUESTION 19. At our Pack forest site, what is the length of 0.005 degrees in meters for longitude and latitude?**
A spatial resolution of hundreds of meters doesn’t sound that good, but considering that GEDI data is aiming for near global coverage, the resolution becomes more reasonable.
Lets plot our height metrics:
?levelplot
plot(rh100metrics$min, main = "rh100 Minimum", xlab = "Longitude (Degrees)", ylab = "Latitude (Degrees)")
plot(rh100metrics$max, main = "rh100 Maximum", xlab = "Longitude (Degrees)", ylab = "Latitude (Degrees)")
plot(rh100metrics$mean, main = "rh100 Mean", xlab = "Longitude (Degrees)", ylab = "Latitude (Degrees)")
plot(rh100metrics$sd, main = "rh100 Standard Deviation", xlab = "Longitude (Degrees)", ylab = "Latitude (Degrees)")
Now we can save these as rasters for viewing in ArcGIS Pro
crs(rh100metrics) <- "EPSG:4326"
writeRaster(rh100metrics$min, "Lab9/Lab9Data/rh100_min.tif", overwrite = TRUE)
writeRaster(rh100metrics$max, "Lab9/Lab9Data/rh100_max.tif", overwrite = TRUE)
writeRaster(rh100metrics$mean, "Lab9/Lab9Data/rh100_mean.tif", overwrite = TRUE)
writeRaster(rh100metrics$sd, "Lab9/Lab9Data/rh100_sd.tif", overwrite = TRUE)QUESTION 20: Submit 4 screenshots of your four rh100 plots. Caption them with a full description of what they represent.
Compare to lidar…
From the Washington DNR lidar portal, grab the DTM from PackForest if you don’t already have it. Just the one file, should be 122.87MB.
Bring the DTM into R:
PFDTM <- rast("Lab9/Lab9Data/pack_forest_2013_dtm_1.tif")Lets check how well the lidar DTM agrees with the elevation values that were captured in the GEDI level 1b data. We already have the level 1b data in a spatial data frame so we can overlay the points across the lidar DTM and extract the DTM values at each point.
This code extracts the DTM values with the GEDI datapoints we
converted to an sf spatial object
level1BGeo_sf_prj <- st_transform(level1BGeo_sf, crs(PFDTM))
GEDI_DTM <- terra::extract(PFDTM, level1BGeo_sf_prj, bind=TRUE)
GEDI_DTM <- na.omit(GEDI_DTM) #remove all NA rowsNow we can plot the elevation from GEDI and the elevation from the DTM together
plot(GEDI_DTM$elevation_bin0, GEDI_DTM$pack_forest_2013_dtm_1)
Let’s evaluate further and make a linear model
reg <- lm(GEDI_DTM$pack_forest_2013_dtm_1 ~ GEDI_DTM$elevation_bin0)
summary(reg)##
## Call:
## lm(formula = GEDI_DTM$pack_forest_2013_dtm_1 ~ GEDI_DTM$elevation_bin0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -432.83 -87.02 -10.53 89.40 428.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.5423 22.0984 2.378 0.0178 *
## GEDI_DTM$elevation_bin0 2.8663 0.0524 54.699 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 143.4 on 557 degrees of freedom
## (1621 observations deleted due to missingness)
## Multiple R-squared: 0.8431, Adjusted R-squared: 0.8428
## F-statistic: 2992 on 1 and 557 DF, p-value: < 2.2e-16Now let’s make a new plot with the regression line
plot(GEDI_DTM$elevation_bin0, GEDI_DTM$pack_forest_2013_dtm_1, xlab = "GEDI Elevation", ylab="ALS Elevation")
abline(reg, lty=3, col ="red", lwd = 3)
abline(0, 1, col = "blue", lty = 5, lwd=2)
QUESTION 22: Submit a screenshot of the scatterplot but using different lty and lwd. Report adjusted R2. Provide a full caption of the figure.
LAST STEP! With your ALS DTM and your GEDI_1B_GEO point layers in ArcGIS, change the colors of the GEDI_1B_GEO symbols to be graduated colors with the elvtn_0 field.
QUESTION 23: Create a map with scale bar and north arrow of pack_forest_2013_dtm and the GEDI_1B_GEO points. Fully caption the figure.